home *** CD-ROM | disk | FTP | other *** search
/ Our Solar System / Our Solar System.iso / miscprog / librat / epsiln.c < prev    next >
C/C++ Source or Header  |  1990-06-20  |  2KB  |  66 lines
  1. /* Obliquity of the ecliptic at Julian date J
  2.  *
  3.  * IAU Coefficients are from:
  4.  * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
  5.  * "Expressions for the Precession Quantities Based upon the IAU
  6.  * (1976) System of Astronomical Constants,"  Astronomy and Astrophysics
  7.  * 58, 1-16 (1977).
  8.  *
  9.  * Before or after 200 years from J2000, the formula used is from:
  10.  * J. Laskar, "Secular terms of classical planetary theories
  11.  * using the results of general theory," Astronomy and Astrophysics
  12.  * 157, 59070 (1986).
  13.  *
  14.  *  See precess.c and page B18 of the Astronomical Almanac.
  15.  * 
  16.  */
  17.  
  18. #include "prec.h"
  19. /* The results of the program are returned in these
  20.  * global variables:
  21.  */
  22. DOUBLE jdeps = -1.0; /* Date for which obliquity was last computed */
  23. DOUBLE eps = 0.0; /* The computed obliquity in radians */
  24. DOUBLE coseps = 0.0; /* Cosine of the obliquity */
  25. DOUBLE sineps = 0.0; /* Sine of the obliquity */
  26. extern DOUBLE eps, coseps, sineps, STR;
  27. DOUBLE SIN(), COS(), FABS();
  28.  
  29. void epsiln(J)
  30. DOUBLE J; /* Julian date input */
  31. {
  32. DOUBLE T;
  33.  
  34. if( J == jdeps )
  35.     return;
  36.  
  37. T = (J - 2451545.0)/36525.0;
  38.  
  39. /* This expansion is from the AA.
  40.  * Note the official 1976 IAU number is 23d 26' 21.448", but
  41.  * the JPL numerical integration found 21.4119".
  42.  */
  43. if( FABS(T) < ((DOUBLE )2.0) )
  44.     eps = (((1.813e-3*T - 5.9e-4)*T - 46.8150)*T + 84381.448)*STR;
  45.  
  46. /* This expansion is from Laskar, cited above.
  47.  * Bretagnon and Simon say, in Planetary Programs and Tables, that it
  48.  * is accurate to 0.1" over a span of 6000 years. Laskar estimates the
  49.  * precision to be 0.01" after 1000 years and a few seconds of arc
  50.  * after 10000 years.
  51.  */
  52. else
  53.     {
  54.     T /= 10.0;
  55.     eps = ((((((((( 2.45e-10*T + 5.79e-9)*T + 2.787e-7)*T
  56.         + 7.12e-7)*T - 3.905e-5)*T - 2.4967e-3)*T
  57.     - 5.138e-3)*T + 1.99925)*T - 0.0155)*T - 468.093)*T
  58.     + 84381.448;
  59.     eps *= STR;
  60.     }
  61.  
  62. coseps = COS( eps );
  63. sineps = SIN( eps );
  64. jdeps = J;
  65. }
  66.